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ABSTRACT 


Econometricians must choose between many methods for estimating p, the 
autocorrelation coefficient, in a first order autoregressive process. This thesis examines 
the performance of four estimators in a Monte Carlo simulation. The methods 
examined are Durbin-Watson, Beach-MacKinnon, Theil-Nagar and Prais-Winsten. 
The autocorrelation coeficient, p, was varied from .2 to .9 and each method provided 
estimates of p and β, the regression coefficient, for 1000 replications. The results 
presented here are similar to those found in previous comparisons. Specifically, 
Ordinary Least Squares was found to be an efficient estimator of p when 
autocorrelation is present only to a slight degree. Of the four estimators examined, the 
performance of Theil-Nagar proved superior in estimating both p and p for small 
values of the correlation coeficient. Beach-MacKinnon, on the other hand, while 
containing a large bias in the estimation of f, is the more efficient estimator of D for 


large values of p. 
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I INTRODUCTION 


A. BACKGROUND 

Autocorrelation exists in a regression model when the error terms are no longer 
independent but are correlated. In the examination of time series data autocorrelation 
is a common phenomenon and can lead to problems if Ordinary Least Squares (OLS) 
estimation procedures are used. The purpose of this thesis is to examine and compare 
four different estimates of the autocorrelation coefficient, p, the estimation of which is 
essential to the resolution of OLS deficiencies. The four estimators to be examined are 


the Durbin- Watson, Theil-Nagar, Deach-MacKinnon, and Prais- Winsten. 


D. PROBLEM STATEMENT 

[n the standard regression modcl v^ X[ - e, y is a Tx1 vector of observations of 
a dependent variable, X is a TxK design matrix and f is a Kxl1 vector. The variable e 
15 а Txl vector of unobservable random errors with E(c)=0 and covariance matrix, 
Е(ес)- 6414. Thus, in the standard model, the random vector e contains elements 
which are pairwise uncorrelated with identical means and variances. In the presence of 
autocorrelation this strong assumption is violated. That is, the error terms are no 


longer independent but are correlated. The regression model becomes, 


y XP ин 21:25:51 (eqn ΣΡ 
where e, — pe, | t vy, 
"ποπ 
E(vv’) = σ΄]. 


This is known as a first order autoregressive or AR(1) process. As illustrated by 
equation l.l, e, 1s expressed linearly in terms of the e, and another random error 
term v... The assumption of zero mean and constant variance provides v, with all the 
nice properties of e, in the standard model. This process may occur for a variety of 
reasons, some of which are: 

1 Omitted explanatory variables.. If a correlated, explanatory variable has. been 
excluded from the, design matrix its exclusion will be reflected in the correlation 
of the random variable с. 


2 20 ication of the mathematical form of the model. \f the wrong mathematical 
relationship is chosen the values of e may be dependent. 


3 Interpolations in the statistical observations. If the observational data is smoothed 
autocorrelation may result. 


4 Mispecification of the true random error. Dependence among the error terms may 
occur naturally. [Ref. 1:р. 204] 


Utilizing OLS to estimate the regression coefficient, D, in the presence of an AR(I) 
process can lead to problems. Generally, there are two consequences to consider. The 
first is that the OLS estimator of the coefficients will be unbiased but will not be very 
efficient. The second consequence is that the OLS variance estimator is biased. For 


these reasons it is useful to investigate other methods to estimate PB [Ref. 2:p. 439]. 


C. ESTIMATORS 

When p is known, the process is easily accounted for using Generalized Least 
Squares or Weighted Least Squares methods [Ref. 3]. However, the usual situation 1s 
that p is unknown and must be estimated. A number of methods have been proposed 
to estimate p and properly account for OLS deficiencies in estimating P. Chapter 2 will 


develop the four estimators mentioned above and examine the autocorrelation process. 


D. SIMULATION 

Each of the estimators considered here have the same asymptotic properties 
therefore any decision on which one to use must be based on small sample analysis and 
Monte-Carlo evidence. Therefore, a simulation will be created in which the data is 
generated according to guidelines presented in previous studies with equation 1.1 as the 
model. The actual values of p will be varied from .2 to .9. The four estimation 


techniques will then provide estimates of p and D for 1000 replications. 


E. MEASURE OF EFFECTIVENESS 

To provide an indication of which estimator performs best the mean square error 
of both P and B will be estimated for each estimator. Prior results for different sets of 
estimators indicate that no one estimator will prove superior over the entire range of p 


but that one or two may out perform the others over specific intervals. 


II. ESTIMATION 


А. GENERAL 

This chapter attempts to explore the theory behind both the first order process 
and four estimators developed to properly account for it. Three of these 
(Durbin-Watson, Thiel-Nagar, and Prais-Winsten) are categorized as estimated 
generalized least squares estimators. The fourth (Beach-MacKinnon) is a maximum 


likelihood estimator. 


B. PROCESS 


The first order process can be written as 


уре о (eqn 2.1) 
уеге ер к= рес 

БУ) 

E(vv’) = 071, 

Е(у,2)= o, ^, and 

E(v,v.)70 for szt. 


The parameter p is generally unknown and along with f must be estimated. The 
statistical properties of the random error, v, listed in equation 2.1 are identical to those 
listed for e in the general linear model. The statistical properties of e under these new 


assumptions are quite different. Judge [Ref. 4:p. 438] shows that 


Ес... (eqn 2.2) 
(2) 
and 
E(e,^) » 6, ^ — c, ^/(1-p)^ . (сап 2.3) 


The covariance between errors s periods apart is no longer zero and is given by 


Ё(еүсү 27 (е, + ҳе) = (7 (1-95). (eqn 2.4) 


The covariance matrix for e is now easily written as 


Ф = E(cc)= | (eqn ДОУ) 
az | Т-! 
p 1 Л 
η; 
р1:1 1-2, , , | 


or utilizing the following convention, 


@ = AJ (eqn 2.6) 


where M! = 


| p р? ο. 
ρ 1 р2 pl 
р? р 1 қ ; І pi? 
= 10-р?) 
pl pl? pl^ | | | ] | 


Thus, the assumptions made about the error term, e, in the standard linear model 
no longer hold for the autoregressive case. Specifically, due to autocorrelation the error 


covariance matrix 1s no longer written as c^] but is now оу 


When an attempt is made to perform a least squares fit to the data in the presence 
of an AR(1) process there are two problems to consider. 


: s | | 
1 ын least squares estimator = (X XYİX'y will be unbiased but will not be very 
efficient. 


° ° А ”хлүэ 5 A 7 , 7 J£" л 
zaslhe least squares covariance matrix G^(X > l with F = (y-Xb) (y-Xb)/(T-K) 
will be a biased estimator of the variance of p. 


In the presence of positive autocorrclation Judge [Ref. 4:p. 439] notes that with 


OLS estimation the bias of the standard error of p will very likely appear as an 
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underestimate. Park and Mitchell [Ref. 5:p. 16] warn that О ο. 
underestimates the variance of D for p > 0.4. This understatement makes the 
estimates themselves appear much more significant than they actually are and makes 


hypothesis testing of the slope coefficients unreliable. 


C. METHODS OF ESTIMATION 
l. Generalized Least Squares Estimation 
When apriori information is available about ‘F, the most convenient estimate 
for the regression coefficient, фр, is obtained by applying least squares estimation 


techniques to the transformed model, 


ее" (eqn 2.7) 
where Y” = py 
A Uum PX 
cape 


The transformation matrix P is the TxT matrix 


рО 0 0 
-p | ο 0 0 
0 -p 1 0 0 
Р = 0 ПЕ” l 0 
0 0 0 -p ] 
where PP = HL. 


This method is known as the Generalized Least Squares (GLS) estimation. 
2. Estimated Generalized Least Squares 
The usual case is that p is unknown and must be estimated. Once an estimate 
for p (р) Is computed one can substitute p into the P matrix and proceed with the GLS 
method outlined above. This is known as Estimated Generalized Least Squares 
(EGLS) estimation. The computational form of the alternative estimators for p 


discussed are as follows: 


a. Durbin-lVatson 
The statistic 


-y 
4 


Sh 
d= X (e,-0.,)7/),8,  t=1,.,T (eqn 2.8) 
172 ~ των 
уућеге е, = урХ 


СИ vised. to test for first order autoregressive errors. As the number of 
observations (T) increases it can be demonstrated that d approaches the least squares 


estimator of f or 


= 1- (9/2). (egn 2.9) 


The Durbin-Watson statistic is provided by most least squares computer packages and 
1s very easy to use. It also is an example of a two-stage estimator. That is, it first 
estimates the correlation parameter and then uses this estimate to compute the 
generalized least squares estimates for f. 
b. Theil-Nagar 
A modification of the Durbin-Watson estimator suggested by Ilenri 
Theil and A. L. Nagar is 


A : a Р 
p= (T(1-(d/2)) + K?) Т2 - К^). (eqn 2.10) 
Theil and Nagar claim that this estimator is an improvement over Durbin-Watson if 


the first and second differences of the explanatory variables are small when compared 


to their corresponding ranges [Ref. 6]. Like Durbin-Watson, it also is a two-stage 


estimator. 
c. Prais-Winsten 
A minimum sum of squares approach to estimating p yields, 
A T xt T- X | 
P=) &(0,1/2, $2. t7 1,..,1 (eqn 2.11) 
= = 
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This estimator can be employed in both a two step and an iterative procedure. This 
paper, however, considers only the following iterative form: 
l. бег = 0. 
2. Transform the variables in accordance with the transformation matrix and 
equation 2.7. A 
3. Calculate the least squares estimate of | ST orar on p. 
4. Calculate the estimate of p conditional on B by using equation 2.11. 
5. If the absolute difference in f from the previous iteration 1s sufficiently small (less 
than 0.00001) stop. If not goto step zZ) Re pam 


3. Maximum Likelihood Estimation 
A maximum likelihood (ML) estimator is the value of O which maximizes the 
value of the likelihood function L(0). Under the assumption that Y has a multivariate 
normal distribution with mean Xp and covariance matrix o Y, the likelihood function 


IS 


г(В,р,с2) = С- (1/20 Ху-Х УУР” Цу-Х) (eqn 2.12) 
where C = -(T/2)Inc^, -- (1/2)In (1-02) . 


2 


The ML estimators for D, p, and c, ^ are those values for which, 


д1/Әф-0, д1/0р=0, д1/06^,=0. (eqn 2.13) 


Solutions to equations 2.13 are very difficult to derive. Beach and MacKinnon 


2 


[Ref. 8:p. 54] use an ML estimator for 6,“ and substitute into equation 2.12. The 


result is the concentrated likelihood function, 
L(B.p) » K-(T/2)n((y-XD) 1H y-Xpy(1-p2)!/ 1) (eqn 2.14) 
where K 7 (T/2)In( T)-(T/2) . 


They suggest maximizing L(D,p) with respect to p with p held constant and then to 
maximize with respect to p with f held constant. An algorithm to derive this ML 


estimate 1s 
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|. Set p = 0. 

2. Transform the variables in accordance with equation 2.7. 

3. Calculate the least squares estimate of D conditional on р: 

4. Calculate the ML estimate of p conditional on f by solving a cubic equation of 
the untransformed residuals. (see [Ref. 8] for details) 

5. If the absolute difference шр from the previous iteration is sufficiently small (less 
than 0.00001) stop. If not, go to step 2 [Ref. 7]. (iNote: The same procedure was 
employed for iterative Prais- Winsten method except that equation 2.11 was used to 


estimate p.) 


This is not a comprehensive listing of all available estimators for a first order 


process. Other estimators are listed in Judge [Ret. 4]. 


lll. COMPARISON 


A. GENERAL 

The finite sampling properties of the estimators listed here have not been derived. 
Choice of which estimator to use might be based on evidence obtained from Monte 
Carlo simulations. This chapter explains a simulation used and provides a synopsis of 


comparisons reported in the literature. 


B. PREVIOUS COMPARISONS 

There have been a number of studies of estimators for p . Each has concluded 
that OLS has serious deficiencies in the presence of autocorrelation. The majority of 
these papers have settled on two points. First, particularly in small sample sizes 
(T<50) it is best to use estimators that consider all T observations Rao and 
Grilitches concluded that using estimators such as Cochrane-Orcutt that ignore the 
first observation can lead to a substantial loss of efficiency [Ref. 9:p. 269]. These 
results were further substantiated by Beach and Mackinnon. In an attempt to develop 
a computationally efficient algorithm to maximize the likelihood function they 
discovered (for p= 0.6, 0.8, 0.99) significant gains in efficiency to be made by 
employing the first observation. Some of these gains are in the neighborhood of 700 
percent [Ref. 8:p. 55]. Park and Mitchell concluded that retention of this ІШЕ 
observation substantially reduces the risk of collinearity as p approaches 0.9 [Ref. 5:p. 
10]. Kobayashi verified theoretically the experimental results of Park and Mitchell. By 
computing the asymptotic variances of several estimators he demonstrated that the loss 
of efficiency of the Cochrane-Orcutt method was due primarily to ignoring the first 
observation [pet 1ο, 5 


The second point is that the Prais-Winsten solution techniques outperform many 
comparable estimators of the correlation parameter. Spitzer concluded that 
Prais-Winsten “appeared to be the best of all the two stage estimators.” [Ref. 1 1:p. 
44]. Park and Mitchell in a later study comparing Beach-MacKinnon with the iterative 
Prais-Winsten estimator concluded that the iterative Prais-Winsten performs 


“appreciably better in estimating the autocorrelation coeflicient p” [Ref. 7:р. 5]. 
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Although there were no studies found specifically comparing the four estimators 
presented here, each has demonstrated a superiority to OLS in the presence of a first 


order process. 


C. MODEL AND DATA GENERATION 
Equation 2.1 was utilized as the model with the first term in the vector e 


generated in the following fashion, 
ei У 1-р) 2. (сап 3.1) 


In order to conform with previous comparisons, the data utilized in this 
experiment is identical to that used in Beach and MacKinnon [Ref. 8]. Two sample 
sizes of 20 and 50 observations were used. The untrended explanatory variable, X, was 
p Was drawn from Х(0, 0.0056). 


Although autocorrelation in theory may be positive or negative, in econometric data it 


drawn from N(0, 0.0625) and the random error, v 


is almost always positive [Ref. 1:p. 201]. For this reason p was varied from 0.2 to 0.9. 


D. VALIDATION 

The data generation program was checked to ensure the normality of e using the 
Chi Square Goodness of Fit test. The normality assumption was accepted at a 0.3684 
level. Finally, in order to ensure each estimator performed properly the random 
portion of the model, specifically the random variable V, was removed. This allowed 
the estimators to function in a deterministic fashion. Data were then generated and 
submitted to each estimator for values of p equal 0.2, 0.6, 0.8. The results are 


presented in Table I , illustrate that the estimators are functioning properly. 


E. SIMULATION 

For each run the values of the regression coefficients, By and PB), were set to | 
and 1. The variables X, and v, were drawn from the normal distributions discussed 
earher. The dependent variable y, was calculated using equation 2.1. Since the 
ultimate objective was to generate residuals to send to the four estimation routines, a 
regression was then performed of y on X and residuals calculated using, 


^ 


E Е В 
ер = у, - хВ аа №. (eqn 3.2) 


TABLE 1 
ESTIMATES OF AUTOCORRELANONCOBERIEON 





The values of the residuals were then sent to each estimation routine. Estimates of 
[| and p were determined for values of p equal to .2, .3, .4, .5, .6, .7, .8, and .9. Bach 


estimate was replicated 1000 times for the sample sizes of 20 and 50. 


F. NIEASURES OF EFFECTIVENESS 
In order to compare the performances of the estimators, two MOE’s were used. 
The mean square error (MSI) of p was estimated for each estimator. This represents 
the expected squared error made in estimating p. The following computational form of 
MSE was used, 
"DAL 


Y (p-p)2/1000  i=1,...,1000. (eqn 3.3) 


/5/ 


| ^. 
The successive values of MSE of p were then plotted against the actual p to 


examine performance over the range of p. 


The second MOE examined the relative efficiencies of the regression coefficient a 
defined in [Ref. 7:p. 7]. A ratio of MSE of [ for a particular estimate to the MSE of B 


for the OLS estimate allows the examination of the relative gains in using particular 
techniques over OLS. Since the proper estimation of f is paramount the efficiency of B 


is predetermined to be the most important MOE. 


IV. RESULTS AND CONCLUSIONS 


A. GENERAL 

The major emphasis of this thesis was to examine the performance of four 
estimators of the autocorrelation coefficient, p, for a first order autoregressive process. 
The estimators examined were Durbin-Watson, Theil-Nagar, Prais-Winsten, and 
Beach-MacKinnon. 

A Monte-Carlo simulation was performed for the following values of p: .2, .3, .4, 
.5, .6, .7, .8, апа .9. Each run was replicated 1000 times for sample sizes of 20 and SB 
The results are recorded in Table IHI. Irrespective of sample size, each of the methods 
underestimate the true value of p but as the number of observations rs increased from 
20 to 50 the bias reduces. As was expected, no one estimator uniformly outperforms 
the others. In both sample sizes, the two stage estimators (Durbin- Watson and 
Theil-Nagar) achieve better results for small p. As the value of p increases, the 
iterative methods (Prais-Winsten and Beach-Mackinnon) perform best. With T=20 


this transition occurs at p=.6 while at 50 observations it occurs earlier at p=.4. 


The discussion of the results will be divided into two sections. The measures of 
effectiveness, as defined in Chapter 3 will first be apphed to the simulation results for 
T=20. This will be followed by an identical approach when the sample size 15 


increased to 50. 


Bae „SANE LESIZEŻ0 
Since performance of an estimator 15 roughly indicated by its mean and variance, 

mean square error (MSE) of each э over the entire sample size was estimated. [he 

results of these calculations are presented in Table It] along with a plot of Nise of ^ 
versus actual values of p in Figure 4.1 . They again indicate that the Theil-Nagar and 

Durbin-Watson estimators are better for smaller values of p (p «.6) and as p increases 

the Prais-Winsten p emerges as the best. On the basis of Figure 4.1 alone, 

Beach-MacKinnon's performance is clearly inferior. However, ш examining the 

efliciency of each estimator in Table IV, Beach-MacKinnon proves to be the most 

efficient in estimating D over the widest range of p. The tie in Figure 4.1 between 


Theil-Nagar and Durbin-Watson 1s resolved in Table IV with Theil-Nagar proving to 


TABLE II 
ESTIMATES OF AUTOCORRELATION COEFFICIENT 


Sample Size 20 


р DW TN PW BM 
2 158 162 113 107 
5 234 239 205 193 
4 310 316 296 279 
5 385 392 387 365 
6 460 467 478 450 
J 533 542 567 515 
8 603 .613 655 617 
9 667 681 741 697 
Sample Size 50 
р DW TN РУУ BM 
2 NIS 161 170 178 
3 279 05 270 270 
4 260 340 360 359 
5 451 432 463 453 
6 544 523 559 547 
7 630 610 651 640 
8 ‚726 ‚700 747 731 
9 812 ‚796 839 820 
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be uniformly more efficient than Durbin-Watson. Table IV also demonstrates that for 


p =.2 OLS 1s at least as efficient as three of the four estimators. 


C. SAMPLE SIZE 50 

The results of the MSE calculations for T=50 are recorded in Table V along 
with a plot of MSE of p versus the actual values of p in Figure 42 ANIĘ 
Durbin- Watson and Theil-Nagar estimators again perform the best for smaller values 
of p (p «.4) and as p increases the Beach-McKinnon and Prais-Winsten estimators of 


p contain the smallest MSE. 


Once again even though the Prais-Winsten p has a smaller MSE than 
Beach-MacKinnon, Table VI illustrates that Beach-McKinnon is a uniformly more 
efficient estimator of the slope coefficient. For the smaller values of p (PSO 
Theil-Nagar 1s more efficient than Durbin-Watson. Table VI also illustrates that OLS 


is at least as efficient as any of the other estimators when p is small. 


D. SUMMARY 

As was found in previous studies when autocorrelation is present only to a slight 
degree (p<.2) the OLS estimator provides an efficient estimate for the regression 
coefficient, D. As the process becomes more significant however, all the estimators 
outperform the OLS solution. In both sample sizes the performance of Theil-Nagar 
and Durbin- Watson are nearly identical with respect to the MSE of p. However, when 
efficiency of the slope coefficient estimate is examined, Theil-Nagar proves to be the 
better 2 stage estimator. Park and Mitchell [Ref. 7:p. 4] found that Prais-Winsten 
performs better in estimating |. The results presented here tend to dispute that 
finding. For while Prais-Winsten has a uniformly smaller MSE of р, 
Beach-Mackinnon provides the most efficient estimator of B . Spitzer, on the other 
hand [Ref. Il:p. 44], which ranked two stage estimators as being the best for values of 
p between .2 and .5, mirrors the results produced here. Apriori knowledge of the 
neighborhood of p will be helpful in selecting the appropriate estimation method. For 
both sample sizes Theil-Nagar appears to be the best for small values of p. 
Beach-MacKinnon, while containing a larger bias for p than does Prais-Winsten, is a 


much more efficient estimator of the slope coefficient for larger values of p. 
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ligure 4.11 Estimated mean square error of p vs. p (sample size = 20). 


TABLE III 
DATA PRESENTED IN FIGURE 4.1 


Sumple Size 20 

MSEDW MSETN MSEPW MSEBM 
1404 1485 боо] 8594 
.6268 .6244 ‚1004 ‚1115 
‚5156 2124 .5621 .5787 
.4159 .4126 .4400 .4604 
2218 „3250 3342 ‚3566 
οι .2495 .2433 .2662 
.1891 .1860 .1698 .1916 
.1407 SN 11141 „158 


p 
22 
К 
4 
3 
.6 
1] 
2 
9 
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Τμ 
EFFICIENCY OF REGRESSION СОГОТ МИСУ ИИИ 


Sample Size 20 
MSE (DW) MSE (TN) MSE[ (PW) MSEB (BM) 


MSEB (OLS)  MSEB (OLS)  MSED(OLS) = MSEB (OLS) 
1.004 9794 1.035 1.041 
.9228 8967 9442 9515 
8218 ‚7929 8325 8342 
7082 6751 7024 6959 
5864 5484 „5652 5515 
4610 4207 4329 4135 
3359 3020 3093 2970 
2251 2087 2077 1892 
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Р ^ ; 
Figure 4.2. Estimated mean square error of f vs. p (sample size 


рТ БЭР УЛ ЭЛ М1610717К12-4-2 


Sample Size 50 


p 
2 
3 
4 
25 
6 
; 
9 
9 


MSEDW MSETN 


‚1010 
‚3500 
4383 
192088 
2335 
.1578 
.0906 
.0500 


TABLE V 


.6766 
530 
4106 
ДЕЙ 
2220 
.1526 
.0942 
.0509 
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MSEPW MSEBM 
‚1065 
2578 
3315 
3156 
2215 
‚1449 
0851 
0417 


1055 
2055 
4206 
.3056 
.2116 
215357 
0773 
0360 





50). 





TABLE VI 
EFFICIENCY OF REGRESSION COEFICIENTES TEA 


Sample Size 50 
MSEB (DW) MSEB (TN) MSEB (PW) MSEB (BM) 


MSEB (OLS)  MSEB(OLS) MSEB (OLS) MSEB (OLS) 
1.073 1.041 1.016 1.058 
9985 9482 9714 9562 
8850 8255 8635 8255 
1452 6859 посо 6825 
5920 5420 5870 5406 
4366 4020 4453 4020 
‚2889 2690 3067 2700 
1589 1505 1738 1505 
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APPENDIX 
PROGRAM LISTINGS 


This appendix contains listings of the programs utilized in the analysis performed 
herein. All of the functions are written in FORTRAN and contain the necessary 
documentation. The Monte Carlo simulation was performed using the Advanced 
Simulation and Statistics Package [Ref. I2] developed by P. A. Lewis. Since the 
package only allows for the simultaneous comparision of 3 estimators, 2 functions were 
Mevcloped Гог each sample size. The first, SIMS generates estimates for 
Durbin-Watson, Theil-Nagar, and Prais-Winsten for a sample size of 20. SIMSA 
meanwhile, generates estimates for Beach-MacKinnon for the identical sample size. 
Routines for Durbin-Watson and Theil-Nagar were included in SIMSA to ensure the 
results were comparable to SIMS. SIMSB and SIMSC perform іп a simular fashion for 
sample size of 50 and therefore were not included. The Advanced Simulation and 
Statistics Package computes the mean square error ог Ф for each estimator 
automatically. The mean square error for the b estimates was accomplished by the 
MSEB function. 


SIMS 
DIMENSION EHAT(20) 
COMMON /MYDATA/ К,Т,АМ5,Ү1,Х 
COMMON /DATA1/ IX1A,RHO 
REAL*4 Y(5000) , YMIN, YMAX, PMEAN( 3) 
CHARACTER*80 T1,T2,T3 
INTEGER N,M,NE(8),L,D,RG,SEI,SVS,NEST,NSR, IXI, IX2, IX3 
EXTERNAL DATGEN, DURWAT, BEAMAC, PRAWIN, LSEB, DCALC, TRANSF 
EXTERNAL LNORM,SIMTBD,GMPRD 
. NRz20 
T=20 
K=2 


10 


120 


115 


61 


999 


OPEN (UNIT=19, FILE='MONICA') 
OPEN (UNIT=21, FILE='MARGE! ) 
OPEN (UNIT=51,FILE='AMBROSE") 
OPEN (UNIT=41,FILE='DAT2') 

OPEN (UNIT=61,FILE='DAT3') 

READ (19,*) ANS 
READ(19,*,END=999) N,M,L,D,RG,SEI,SVS,NEST,NSR 

READ(19,*)YMIN, YMAX 

READ(19,*) (NE(I),I=1,L) 

READ(19,120) IX1,IX2,IX3 
FORMAT(I5,1X,I5,1X, I5) 

READ (19,115) T1 
FORMAT(A80) 

READ(19,115) T2 

READ (19,115)T3 

READ(19,*) (РМЕАМ(1),1-1,3) 

READ(19,*) RHO 

READ(19,61)IX1A 
FORMAT(15) 


CALL FOR SIMTBD 


CALL SIMTBD (IX1,IX2,IX3,Y,N,M,NE,L,D,NSR,RG,SEI,SVS, 
*YMIN, YMAX , NEST , DATGEN , DURWAT , T 1, DATGEN , BEAMAC , T2, DATGEN , PRAWIN,T3, 
*PMEAN) 

GO TO 10 
WRITE(6,*)'END OF DATA INPUT! 

STOP 

END 


ΧΧΧ Ae oe oe oe eoe oe * ἈΚ * eoe oe oe coe ve *x * e eee oe e ec o oe oo oe oe oe oe o Ok KKK KKK KK KK KK KK KK KK KKK KK 


жхххжжххжхжхх DATA GENERATION SUBROUTINE XX xo ove eo oo Ἃς * * * ve o e e x ἈΚ Χ 


ΧἀχΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧ 


SUBROUTINE DATGEN (IX1,EHAT,NR) 
DIMENSION BHAT(2),YSTAR( 20) ,R2(20) ,U(20), 
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Co. Co = Ga CJ 


O O 


38 


O O O O 
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> 


*E(20) ,YHAT(20) ,EHAT(20) ,XSTAR(20,2) 
* Y1(20) ,X(20,2),V(20) 
COMMON /MYDATA/ К,Т,АМ5,Ү1,Х 
COMMON /DATA1/ IX1A,RHO 
INTEGER IX1,IX1A,NR 


GENERATE THE RANDOM ERROR 


CALI SNOR (1X1 (USNR, 1,0) 


ADJUST THE VARIANCE OF R. E. IAW BEACH AND MACKINNON( 1978) 
51-11 
V(I)sU(I)*. 06 
CONTINUE 


GENERATE THE ERROR FOR THE STAND LINEAR MODEL 


E(1)=V(1)/(1-(RHO**2))**0. 5 
DO 31 J=2,T 
E(J)=RHO*E(J-1)+V(J) 
CONTINUE 


GENERATE THE EXPLANATORY VARIABLES IAW RAO AND GRILITCHES (1969) 


DO 32 1-1,20 
X(1,1)=1 
CONTINUE 
CHANGE IX1 IN ORDER TO AVIOD COLLINEARITY 
IX1A=IX1+19 
CALL SNOR(IX1A,R2,NR,1,0) 
DO 33 J=1,20 
X(J,2)=R2(J)*. 25 
CONTINUE 
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0 С> OQO CO ΕΘ. τμ 


5ο 


100 


СЭ сэ сэс 


50 


THE TRUE BETA ОВАА 


GENERATE THE INDEPENDENT VARIABLE 


DO 35 I=1,20 
Y1(T)=(X(T,1)+X(1,2))+EC(1) 
CONTINUE 


GENERATE THE LEAST SQUARES ESTIMATOR FOR BETA 


CALL LSEB(X,Y1,BHAT) 
PRINT LSEB TO A FILE 

IF(ANS .EQ. 2) WRITE(61,201)BHAT 
FORMAT(F11. 8,2X,F11. 8) 


GENERATE YHAT 


DO 100 1=1,20 
YHAT(1)=X(I,1)*BHAT(1)+X(1,2)*BHAT(2) 
CONTINUE 


GENERATE EHAT 


DO 50 I=1,20 
ЕНАТ(1)-Ү1(1)-ҮНАТ(1) 
CONTINUE 


RETURN 
END 


Ἂ Χ Χ x x x x x x x x x GG GG X xk *< DURBIN WATSON X Ἂς x x x x x x x x x x x x x x x x x x x x x x x x 
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C 


701 
C 
C 
C 


ἘΠ СЭ?" СЭ. СЭ 


801 


THIS FUNCTION COMPUTES THE DURBIN-WATSON ESTIMATE OF RHO 
REAL FUNCTION DURWAT (EHAT,NR,WI) 
DIMENSION EHAT(20),X(20,2) , Y1(20) ,XSTAR1( 20,2) , YSTAR1(20) , BHAT1(2) 
COMMON /MYDATA/ K,T,ANS, Y1,X 
CALL DCALC (EHAT,T,D) 
DURWAT=1-D/2 


CALL TRANSF(X ,Y1,DURWAT ,XSTAR1,YSTAR1) 

CALL LSEB (XSTAR1,YSTAR1,BHAT1) 

IF (ANS .EQ. 1 ) WRITE(21,701) ВНАТ1 
FORMAT(F11. 8,2X,F11. 8) 


RETURN 
END 


Ἂ KO Ce oe Xe eee eee < < < x *< THEIL NAGAR Ἂς Ἂς Ἂς eo Ἃς Ἃς ἋΚ ἋΚ Ἃς Ἃς ἋΚ Ἃς eoe vede x x * x Ἃς Ἃς Ἃς Ἃ Ἃς Χ 


THIS FUNCTION COMPUTES THE THEIL-NAGAR ESTIMATE OF RHO 
REAL FUNCTION THENAG (EHAT,NR,WI) 
DIMENSION EHAT(20) , YSTAR2(20) ,XSTAR2(20,2) ,BHAT2(2) 
* Y1(20),X(20,2) 
COMMON /MYDATA/ K,T,ANS,Y1,X 
CALL DCALC (EHAT,T,D) 
THENAG=( (T**2)*(1-0/2)+K**2)/(T**2-K**2) 
CALL TRANSF(X,Y1,THENAG, XSTAR2,YSTAR2) 
CALL LSEB (XSTAR2,YSTAR2 ,BHAT2) 
IF (ANS . EQ. 1 ) WRITE(31,801) BHAT2 
πο DUI ο 2X,F11.8) 
RETURN 
END 
*kkkkkkkkkkkkkkKkkkkkkKkkK DRAIS WINSTEN ***x*% &*# @% Y W% xxx 
THIS FUNCTION COMPUTES THE PRAIS-WINSTEN ESTIMATE OF RHO 
REAL FUNCTION PRAWIN(EHAT,NR,WI) 


~ 


21 


25 


83 


DIMENSION EHAT3(20), YHAT3(20) , YSTAR3( 20) ,BHAT3(2) , 
*EHAT(20) ,XSTAR3(20,2) 
* Y1(20),X(20,2) 
COMMON /MYDATA/ К,Т,АМ5,Ү1,Х 
N-0 
RHO3=0 
N=N+1 
CALL TRANSF (X,Y1,RHO3,XSTAR3, YSTAR3) 
CALL LSEB (XSTAR3,YSTAR3, BHAT3) 
GENERATE YHAT3 
DO 83 I=1,20 
YHAT3(1)=X(I ,1)*BHAT3(1)+X(I ,2)*BHAT3(2) 
CONTINUE 
DO 4 I=1,T 
EHAT3(I)=Y1(1)-YHAT3(1) 
CONTINUE 


RHONUM=0 
RHODEN=0 
DO 5 I=2,T 
RHONUM=RHONUM+( EHAT3(1)*EHAT3(1-1)) 
CONTINUE 


DO 6 I=2,T-1 
RHODEN=RHODEN+( EHAT3(1)**2) 
CONTINUE 
PRAWIN=RHONUM/RHODEN 
CHECK FOR PRAWIN WHICH ARE OUT OF BOUNDS 
IF(PRAWIN. GE. 1)THEN 
PRAWIN=0. 99999 
ELSE IF (PRAWIN. LE. -1)THEN 
PRAWIN=-0. 99999 
END IF | 
COMPARISION OF RHO3 AND PRAWIN IF DIFF .LT. 0.0901 THEN END 
IF(ABS( RHO3-PRAWIN). GT. . 0001) THEN 
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C01 


ο ο ον ον το το 


41 
40 


RHO3=PRAWIN 
GO TO 98 
ELSE 
PRAWIN=PRAWIN 
END IF 
IF (ANS .EQ. 1 ) WRITE(41,901) BHAT3 
FORMAT(F11. 8,2X,F11. 8) 
RETURN 
END 


THE FOLLOWING SUBROUTINES AID IN THE COMPUTATION OF THE FOUR 
ESTIMATORS OF RHO. 


X Yo +++... SUBROUTINE | SEB ΧΧΧΧΧΧΧΧΧἈΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧ ΧΧΧ 


ο ΗΕ ESEB WICE COMPUTE THE LSE OF B 


SUBROUTINE LSEB(X,Y1,BHAT) 
DIMENSION BHAT(2) ,Y1(20) ,X(20,2),XTRNSP(2,20),XI(2,2),H(2,20), 

*XPRIX(2,2) 

X TRANSPOSE 
DO 40 I=1,20 

DO 41 J=1,2 
XTRNSP(J,I)=X(I,J) 
CONTINUE 

CONTINUE 

MULTIPLY X TRANSPOSE AND X 
CALL GMPRD(XTRNSP,X,XPRIX,2,20,2) 

CALCULATE INVERSE OF X PRIME X 
DETR=1/(XPRIX(1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 
XI(1,1)=DETR*XPRIX(2,2) 

XI(1,2)=DETR*(-XPRIX(1,2)) 
XI(2,1)=DETR*(-XPRIX(2,1)) 
XI(2,2)=DETR*XPRIX(1,1) 

MULTIPLY INVERSE AND TRANSPOSE 
CALL GMPRD(XI,XTRNSP,H,2,2,20) 

DO 99 I=1,2 


33 


99 


Су С> СУ с С? 


BHAT(I)=H(I,1)*Y1(1)+H(I,2)*Y1(2)+H(I,3)2Y1(3) 
*+H(1,4)*Y1(4)+H(1,5)*Y1(5) 
*+H(1,6)*Y1(6)+H(I,7)*Y1(7)+H(1,8)*Y1(8)+H(1,9)*Y1(9) 
*+H(1,10)*Y1(10)+ 
*H(I,11)*Y1(11)*H(I,12)*Y1(12)4H(I, 13) *Y1(13)4H(I, 14) *Y1(14) 
*+H(I,15)*Y1(15)+ 
*H(1,16)*Y1(16)+H(I,17)*Y1(17)+H(1,18)*Y1(18)+H(1,19)*Y1(19)+ 
*H(I,20)*Y1(20) 

CONTINUE 

RETURN 

END 


KKKKKKKKKKKKK SUBROUTINE NEALE KEKKKKAKK KKK KK KKK KK KKK KK KK KKK 


SUBROUTINE DCALC WILL COMPUTE THE DURBIN STATISTIC D 


SUBROUTINE DCALC(EHAT,T,D) 
DIMENSION 01(20),02(20),ЕНАТ(20) 
DNUM=0 
DDEN=0 
DO 1 I=2,T 
D1(I-1)=(EHAT(I)-EHAT(I-1))**2 
DNUM=DNUM+D1(I-1) 
CONTINUE 
DO 2 J=1,T 
D2(J)=EHAT(J)**2 
DDEN=DDEN+D2(J) 
CONTINUE 
D=DNUM/DDEN 
RETURN 
END 


KKKKKKKKKKKKKKKKK SUBROUTINE TRANSF жжжххххххххххххххххххх 


SUBROUTINE TRANSF IS DESIGNED TO TRANSFORM THE X'S AND Y'S 
ACCORDING TO THE LEAST SQUARES RULE. 
SUBROUTINE TRANSF(X,¥1,RHOHAT , XSTAR, YSTAR) 
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ШІ 


15 
12 


DIMENSION Y1(20),YSTAR(20) ,X( 20,2) ,XSTAR(20,2) 
K=2 
T=20 
Y TRANSFORM 
YSTAR(1)=((1-(RHOHAT**2))**0.5)*Y1(1) 
DO 7 I=2,20 
YSTAR(I)=Y1(1)-(RHOHAT*Y1(1-1)) 
CONTINUE 
X TRANSFORM 
XSTAR(1,1)=(1-(RHOHAT**2))**0.5 . 
DO 9 J=2,K 
XSTAR(1,J)=((1-(RHOHAT**2))**0.5)*X(1,J) 
CONTINUE 
DO 11 L=2,T 
XSTAR(L,1)=1-RHOHAT 
CONTINUE 
DO 12 I=2,T 
00 13 J=2,K 
XSTARCI,J)2X(I,J)-RHOHAT*X(I-1,J) 
CONTINUE 
CONTINUE 
RETURN 
END 


55 


C 
C 


10 


105 


120 


115 


61 


SIMSA 

THE PURPOSE OF THIS PROGRAM IS TO RUN COMPUTE THE FOLLOWING 
ESTIMATORS (DW TN BM) FOR A SAMPLE SIZE OF 20 

DIMENSION EHAT(20) 

COMMON /MYDATA/ K,T,ANS,Y1,X 

COMMON /DATA1/ IX1A,RHO 

REAL*4 Y(5000) , YMIN, YMAX , PMEAN( 3) 

CHARACTER*80 T1,T2,T3 

INTEGER N,M,NE(8),L,D,RG,SEI,SVS,NEST,NSR,IX1, IX2, IX3 

EXTERNAL DATGEN, DURWAT, THENAG, BEAMAC, LSEB, DCALC, TRANSF 

EXTERNAL LNORM,SIMTBD ,GMPRD 

NR=20 

T=20 

K-2 


OPEN (UNIT=19, FILE='MONICA' ) 
OPEN (UNIT=51,FILE='AMBROSE' ) 
READ (19,*) ANS 
READ(19,*,END=999) N,M,L,D,RG,SEI,SVS,NEST,NSR 
READ(19,*)YMIN, YMAX 
READ(19,*) (NE(I),I=1,L) 
WRITE (22,105) (NE(I),I=1,L) 
FORMAT(814) 
READ(19,120) IX1,IX2,IX3 
FORMAT(I5,1X,I5,1X,I5) 
READ (19,115) T1 
FORMAT(A80) 
READ(19,115) T2 
READ (19,115)T3 
READ(19,*) (PMEAN(I),I=1,3) 
READ(19,*) RHO 
READ(19,61)IX1A 
FORMAT( 15) 


36 


999 


Со Сэ сә συ 


Сэ со с> сэ 


3] 


CALL FOR SIMTBD 


CALL SIMTBD (IX1,IX2, IX3, Y, N,M,NE, L,D,NSR,RG, SEI ,SVS, 
*YMIN , YMAX ,NEST, DATGEN , DURWAT , T1, DATGEN, THENAG, T2 , DATGEN, BEAMAC, T3, 
*PMEAN) 

GO TO 10 
WRITE(6,*)'END OF DATA INPUT! 

STOP 

END 


ΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧἈΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧΧ 


kkkkkkkkkkkk* DATA GENERATION SUBROUTINE **xx***kkkkkkkkkkkkkkk 
ОООО УО оники ки екн ананан 
SUBROUTINE DATGEN (IXl,EHAT,NR) 
DIMENSION BHAT(2),YSTAR(20) , R2(20) ,U(20), 
*E(20) , YHAT(20) ,EHAT(20) ,XSTAR(20,2) 
* Y1(20),X(20,2) 
COMMON /МҮРАТА/ К,Т,АМ5,Ү1,Х 
COMMON /DATAl/ IX1A,RHO 
INTEGER. IX1,IX1A,NR 


— — 


GENERATE THE RANDOM ERROR 


CALL SNOR (IX1,U,NR,1,0) 


GENERATE THE ERROR FOR THE STAND LINEAR MODEL 


E(1)-U(1)/(1-(RHO**2))**0. 5 
DO 31 Ј=2,20 
E(J)=RHO*E(J-1)+U(J) 
CONTINUE 
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C GENERATE THE EXPLANATORY VARIABLES IAW RAO AND GRILITCHES (1969) 


DO 32 1-1,20 
Х(1,1)-1 
32 CONTINUE 
C CHANGE IX1 IN ORDER TO AVIOD COLLINEARITY 
C IX1A=1X1+19 
CALL SNOR(IX1A,R2,NR,1,0) 
DO 33 J=1,20 
X(J,2)=R2(J)*. 25 
33 CONTINUE 


THE TRUE BETA EQUALS 1,1 


GENERATE THE INDEPENDENT VARIABLE 


DO 35 I-1,20 
YI(I)z(XCI, 1) *X(1,2) )*E( I) 
P Ol SE 


С 
C 
C GENERATE YHAT 
C 


CALL LSEB(X,Y1,BHAT) 
C ВНАТ(1)=1. 3 
C ВНАТ(2)=1. 1 
00 100 1=1,20 
YHAT(I)=X(I,1)*BHAT(1)+X(1,2) *BHAT(2) 
100 CONTINUE 


GENERATE EHAT 


- GC) €) ©) 


38 


DO 50 1-1,20 
EHAT( I)2Y1( I) - YHAT( I) 
50 ^ CONTINUE 
C 
C 
RETURN 
END 


C Ἂς Xe oe * x x x x x x x x x x Xe x x x x x< DURBIN WATSON Ἂ Χ x x x x x x x X x x x x x x x x x x x x x x x x x< 


REAL FUNCTION DURWAT (EHAT,NR,WI) 

DIMENSION EHAT(20),X(20,2),Y1(20),XSTAR1(20,2),YSTAR1(20),BHAT1(2) 
COMMON /MYDATA/ K,T,ANS,Y1,X 

CALL DCALC (EHAT,T,D) 

DURWAT=1-D/2 

CALL TRANSF(X,Y1,DURWAT,XSTAR1,YSTAR1) 

CALL LSEB (XSTAR1,YSTAR1,BHAT1) 


RETURN 
END 


Ἂ Ἂς x x x x x x x x x x * * Ἃς Ἃς ΧΚ ϑς * *< THEIL NAGAR x o Ἂς Ἃς x x x x x x x x x x x x x x x x x x x x x x x x x *< 


CIE" СОДИС 


REAL FUNCTION THENAG (EHAT,NR,WI) 

DIMENSION EHAT(20) , YSTAR2(20) ,XSTAR2(20,2) ,BHAT2(2) 
* Y1(20),X(20,2) 

COMMON /MYDATA/ К,Т,АМ5,Ү1,Х 

CALL DCALC (EHAT,T,D) 

THENAG=( (T**2)*(1-0/2)+K**2)/(T**2-K**2) 

RETURN 

END 


© 


X x x x x x x x x x x x x x x x x x x x x x * BEACH MACKINNON * * x x x x x x x x x x x x x x x x x x x *x x e Χκ 


29 


98 


83 


21 


12 


105 


REAL FUNCTION BEAMAC(EHAT , NR,WI) 
DIMENSION EHAT4(20) , YHAT4(20) , YSTAR4( 20) , BHAT4(2), 
*¥1(20) ,EHAT(20) ,X(20,2),XSTAR4( 20,2) 
COMMON /MYDATA/ K,T,ANS,Y1,X 
N=0 
RHO4=0 
N=N+1 
CALL TRANSF (X,Y1,RHO04 ,XSTAR4,YSTAR4) 
CALL LSEB (XSTAR4, YSTAR4, BHAT4) 
BHAT 4(1)=1. 0 
ВНАТ4(2)-1.0 
GENERATE YHAT4 
DO 83 I=1,20 
ҮНАТ4(І)-Х(1,1)%ВНАТ4(1)-Х(1,2)“ВНАТ4(2) 
CONTINUE 
DO 4 I=1,T 
ЕНАТ4( І)-Ү1( І)-ҮНАТА(Т) 
CONTINUE 
SUM3=0 
SUM2=0 
SUM1=0 
DO 71 І-2,Т 
50М1-50М1-(ЕНАТ4(1)“ЕНАТ4(1-1)) 
CONTINUE 


DO 72 ес 
SUM2=SUM2+( EHAT4( I-1)**2) 
CONTINUE 


DO 73 I=2,T 
SUM3=SUM3+(EHAT4(1)**2) 
CONTINUE 
DENOM=(T-1)*(SUM2-(EHAT4(1)**2)) 
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901 


Az( -(T-2)*SUM1)/DENOM 


B-(((T-1)*(EHAT4(1)**2))- (T*SUM2) - SUM3) /DENOM 


C=(T*SUM1)/DENOM 


SMALP=B-((A**2)/3) 


SMALQ=C-((A*B)/3)+((2*(A**3))/27) 


THETA=ACOS( (SMALQ*(27**. 5))/(2*SMALP*"((-SMALP)**0.5))) 


BEAMAC IS THE ITERATIVE RHO FOR THIS PROCEEDURE 
BEAMAC=(-2*((-SMALP/3)**0. 5))*COS((THETA/3)+(3.1412/3))-(A/3) 
CHECK FOR BEAMAC WHICH ARE OUT OF BOUNDS 
I F(BEAMAC. GE. 1)THEN 
BEAMAC=0. 99999 
ELSE IF (BEAMAC. LE. -1)THEN 
BEAMAC=-0. 99999 
END IF а 2 
COMPARISION OF RHO4 AND BEAMAC IF DIFF .LT. 0.0001 THEN END 
IF(ABS(RHO4-BEAMAC). GT. . 0001) THEN 
RHO4-BEAMAC 
GO TO 98 
ELSE 
BEAMAC-BEAMAC 


. END IF 


IF (ANS . EQ. 2) WRITE (51,901) BEAMAC 
FORMAT(F15. 11) 

RETURN 

END 


THE FOLLOWING SUBROUTINES AID IN THE COMPUTATION OF THE FOUR 
ESTIMATORS OF RHO. 
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K K X k k k X k k k x x K SUBROUTINE LSEB X k K K x K x x K K k k KX X X x x K k X K k X k X X x X x k K 


SUBROUTINE LSEB WILL COMPUTE THE LSE OF B 


SUBROUTINE LSEB(X,Y1,BHAT) 
DIMENSION BHAT(2),Y1(20),X(20,2) ,XTRNSP(2,20) ,X1(2,2) ,H( 2,20), 

*XPRIX(2,2) | 

X TRANSPOSE 
DO 40 I=1,20 

DO 41 J=1,2 
XTRNSP(J,I)=X(I,J) 
CONTINUE 

CONTINUE 

MULTIPLY X TRANSPOSE AND X 
CALL GMPRD(XTRNSP,X,XPRIX,2,20,2) 

CALCULATE INVERSE OF X PRIME X 

DETR=1/(XPRIX(1,1)*XPRIX(2,2)-XPRIX(1,2)*XPRIX(2,1)) 
XI(1,1)=DETR*XPRIX(2,2) 
XI(1,2)=DETR*(-XPRIX(1,2)) 
XI(2,1)=DETR*(-XPRIX(2,1)) 
XI(2,2)=DETR*XPRIX(1,1) 

MULTIPLY INVERSE AND TRANSPOSE — 

CALL GMPRD(XI,XTRNSP,H,2,2,20) 
DO 99 1-1,2 
BHAT(I)eH(I,1)*Y1(1)*H(I,2)*Y1(2)*H(I, 3) * Y1(3) 

*«H(I,4)*Y1(4)*H(I,5)*Y1(5) 

*-H(I,6)*Y1(6)*H(I,7) *Y1(7)*H( 1,8) *Y1(8)4H( 1,9) *Y1(9) 

*«H(I,10)*Y1(10)* 

*HCI,11)*Y1(11)*H(I,12)*Y1(12)*H(I, 13) *Y1(13) *H( 1, 14) *Y1( 14) 

*«H(I,15)*Y1(15)* 

*H(I,16)*Y1(16)*H(I,17) *Y1(17)*H(I,18) *Y1(18)4H(I, 19) *Y1(19)* 

*H(I,20)*Y1(20) 

CONTINUE 
RETURN 
END 


Ἂ Ἂς K x x K x x x K x < SUBROUTINE пПСАГС o oo oe e ee x x x x x x x x x x x x x x x x x x x x x< 
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C SUBROUTINE DCALC WILL COMPUTE THE DURBIN STATISTIC D 


SUBROUTINE DCALC(EHAT,T,D) 
DIMENSION D1(20),02(20),EHAT(20) 
DNUM=0 
DDEN=0 
DO 1 I=2,T 
D1(I-1)=(EHAT(1)-EHAT(I-1))**2 
DNUM=DNUM+D1(I-1) 
1 CONTINUE 
DO 2 J=1,T 
D2(J)=EHAT(J)**2 
DDEN=DDEN+D2(J) 
2 CONTINUE 
D=DNUM/DDEN 
RETURN 
END 


*X x * x * x  * k x * x * x x < *< SUBROUTINE TRANSF KKKKKKKKKKKKKKKKKKKKKK 


SUBROUTINE TRANSF IS DESIGNED TO TRANSFORM THE X'S AND Y'S 
ACCORDING TO THE LEAST SQUARES RULE. 
SUBROUTINE TRANSF(X,Y1,RHOHAT, XSTAR, YSTAR) 
DIMENSION Y1(20),YSTAR(20) ,X(20,2),XSTAR( 20,2) 
K=2 
T=20 
C Y TRANSFORM 
YSTAR(1)=((1-(RHOHAT**2))**0.5)*Y1(1) 
DO 7 I=2,20 
YSTAR(I)=Y1(1)-(RHOHAT*Y1(1-1)) 
7 CONTINUE 
C X TRANSFORM 
XSTAR(1,1)=(1-(RHOHAT**2))**0. 5 
DO 9 J=2,K 
XSTAR(1,J)=((1-(RHOHAT**2))**0.5)*X(1,J) 


το €) €) (9 C) 
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15 
12 


CONTINUE 
DO 11 L=2,T 
XSTAR(L,1)=1-RHOHAT 
CONTINUE 
DO 12 I=2,T 
DO 13 J=2,K 
XSTAR(I,J)=X(1,J)-RHOHAT*X(I-1,J) 
CONTINUE 
CONTINUE 
RETURN 
END 


44 





MSEB 
С THIS PROGRAM IS DESIGNED TO CALCULATE THE MEAN SQUARE ERROR OF 
C THE BETA VECTOR | 
DIMENSION B1(5000),B2(5000),83(5000),B4(5000),B5(5000) ,B6(5000), 
*B7(5000) ,88( 5000) ,89( 5000) ,810( 5000) ,8Х( 5000) ,ВУ( 5000) 
OPEN (UNIT=21,FILE='DAT1') 
OPEN (UNIT=31,FILE='DAT2') 
OPEN (UNIT=41, FILE='DAT3') 
OPEN (UNIT=51,FILE='DAT4') 
OPEN (UNIT=61,FILE='DATS') 


COUNT=1000 
READ(21,900)(B1(I),B2(I), I=1,1000) 
CALL MSEBET (B1,B2,COUNT , XMSEDW) 
READ(31,900)(83(1),B4(1), I=1,1000) 
CALL MSEBET (B3,B4,COUNT,XMSETN) 
READ(41,900)(B5(I),B6(I), I21,1000) 
CALL MSEBET (B5,B6,COUNT , XMSEPW) 
READ(51,900)(B7(I),B8(I), I21,1000) 
CALL MSEBET (B7,B8,COUNT,XMSEBM) 
READ(61,900)(B9(I),B10(I), I=1,1000) 
CALL MSEBET (B9,B10,COUNT,XMSEOLS) 

900 FORMAT (Е11.8,2Х,Е11.8) 
WRITE(6,*)'MSEDW' 


WRITE(6, *)XMSEDW 
C 

WRITE(6,*)' MSETN' 

WRITE(6,*)XMSETN 
C 

WRITE(6,*)' MSEPW' 

WRITE(6,*)XMSEPW 
C 

WRITE(6,*)'MSEBM' 

WRITE(6, *)XMSEBM 
C 
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WRITE(6,*)'MSELS' 
WRITE(6,*)XMSELS 
STOP 

END 


C *k ye ve s e ye e ὑς e e e ve e e e x SUBROUTINE MSEBET * X x x x x x e x K fe x x e ve oe eoe eoe ve eoe He He He KK 


SUBROUTINE MSEBET(BX,BY,AN,XMSEB) 
DIMENSION BX(5000),BY(5000) , SUM( 5000) 
PLACE=0 
DO 901 I=1,AN 
SUM( I)=((BX(1)-1)*(BY(1)-1))**2 
PLACE=PLACE + SUM(I) 

901 CONTINUE 
XMSEB=PLACE/AN 
RETURN 
END 
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